This script provides a mean of supersaturation time for both species
#load the various libraries
library(lme4)
library(lmerTest)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
library(mmtable2)
library(gt)
library(purrr)
library(stringr)
library(tidyr)
Load the data, assign some of the variables as factors
ek_irrad_data <- read.csv("input_data/mean_supersaturation_by_period.csv")
# assign run as a factor
ek_irrad_data$Run <- as.factor(ek_irrad_data$Run)
#assign temperature as a factor
ek_irrad_data$Temperature <- as.factor(ek_irrad_data$Temp...C.)
#assigns treatment as characters from integers then to factors
ek_irrad_data$Treatment <- as.factor(as.character(ek_irrad_data$Treatment))
# make a new column for delta Ek
ek_irrad_data$deltaEk <- (((ek_irrad_data$day9_ek - ek_irrad_data$ek.1)/ ek_irrad_data$ek.1) * 100)
ULVA ____________________________________________________________
Subset the data by species and remove treatment 2.5 because this had problems with tissue sloughing at/near full moon
ulva <- subset(ek_irrad_data, Species == "ul" & Treatment != 2.5)
ulva$treatment_graph[ulva$Treatment == 0] <- "1) 35ppt/0.5umol"
ulva$treatment_graph[ulva$Treatment == 1] <- "2) 35ppt/14umol"
ulva$treatment_graph[ulva$Treatment == 2] <- "3) 28ppt/27umol"
ulva$treatment_graph[ulva$Treatment == 3] <- "5) 18ppt/53umol"
ulva$treatment_graph[ulva$Treatment == 4] <- "6) 11ppt/80umol"
ulva$treatment_graph[ulva$Treatment == 2.5] <- "4) 28ppt/53umol"
Make a histogram of the data for ulva
hist(ulva$supersat_avg, main = paste("Ulva lactuca"), col = "olivedrab3", labels = TRUE)
ulva %>% ggplot(aes(supersat_avg)) +
geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
theme_bw()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
Run model without interaction between the treatments and temperature - supersat_avg is in minutes Take RLC.Order our of the random effects because causing problems of singularity. R2 is same with or without
hsat_model_ulva <- lmer(formula = supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = ulva)
Check residual plots
hist(resid(hsat_model_ulva))
plot(resid(hsat_model_ulva) ~ fitted(hsat_model_ulva))
qqnorm(resid(hsat_model_ulva))
qqline(resid(hsat_model_ulva))
Check the performance of the model, get R2, summary and print table and plot
performance::check_model(hsat_model_ulva)
## Variable `Component` is not in your data frame :/
r.squaredGLMM(hsat_model_ulva)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
## R2m R2c
## [1,] 0.1963057 0.8597237
summary(hsat_model_ulva)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
## Data: ulva
##
## REML criterion at convergence: 2331.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1705 -0.3833 0.0753 0.4654 2.8339
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant.ID (Intercept) 609.2 24.68
## Run (Intercept) 2724.1 52.19
## Residual 704.8 26.55
## Number of obs: 240, groups: Plant.ID, 96; Run, 7
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 472.758 37.549 5.055 12.591 5.21e-05 ***
## Treatment1 -50.383 44.311 5.002 -1.137 0.307
## Treatment2 -47.642 44.311 5.002 -1.075 0.331
## Treatment3 -80.691 44.311 5.002 -1.821 0.128
## Treatment4 -88.422 44.311 5.002 -1.995 0.103
## Temperature27 7.294 7.839 75.122 0.930 0.355
## Temperature30 4.850 7.839 75.122 0.619 0.538
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmn1 Trtmn2 Trtmn3 Trtmn4 Tmpr27
## Treatment1 -0.835
## Treatment2 -0.835 0.993
## Treatment3 -0.835 0.993 0.993
## Treatment4 -0.835 0.993 0.993 0.993
## Temperatr27 -0.104 0.000 0.000 0.000 0.000
## Temperatr30 -0.104 0.000 0.000 0.000 0.000 0.500
plot(allEffects(hsat_model_ulva))
tab_model(hsat_model_ulva, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
| supersat_avg | ||||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | Statistic | p | df |
| (Intercept) | 472.76 | 37.55 | 398.77 – 546.74 | 12.59 | <0.001 | 230.00 |
| Treatment [1] | -50.38 | 44.31 | -137.69 – 36.93 | -1.14 | 0.257 | 230.00 |
| Treatment [2] | -47.64 | 44.31 | -134.95 – 39.67 | -1.08 | 0.283 | 230.00 |
| Treatment [3] | -80.69 | 44.31 | -168.00 – 6.62 | -1.82 | 0.070 | 230.00 |
| Treatment [4] | -88.42 | 44.31 | -175.73 – -1.11 | -2.00 | 0.047 | 230.00 |
| Temperature [27] | 7.29 | 7.84 | -8.15 – 22.74 | 0.93 | 0.353 | 230.00 |
| Temperature [30] | 4.85 | 7.84 | -10.60 – 20.30 | 0.62 | 0.537 | 230.00 |
| Random Effects | ||||||
| σ2 | 704.81 | |||||
| τ00 Plant.ID | 609.21 | |||||
| τ00 Run | 2724.09 | |||||
| ICC | 0.83 | |||||
| N Run | 7 | |||||
| N Plant.ID | 96 | |||||
| Observations | 240 | |||||
| Marginal R2 / Conditional R2 | 0.196 / 0.860 | |||||
Construct null model to perform likelihood ratio test REML must be FALSE
ulva_hsat_treatment_null <- lmer(formula = supersat_avg ~ Temperature + (1 | Run) + (1 | Plant.ID) +
(1 | RLC.Order), data = ulva, REML = FALSE)
ulva_hsat_model2 <- lmer(formula = supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +
(1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_hsat_treatment_null, ulva_hsat_model2)
## Data: ulva
## Models:
## ulva_hsat_treatment_null: supersat_avg ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_hsat_model2: supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## npar AIC BIC logLik deviance Chisq Df
## ulva_hsat_treatment_null 7 2461.2 2485.6 -1223.6 2447.2
## ulva_hsat_model2 11 2396.4 2434.7 -1187.2 2374.4 72.754 4
## Pr(>Chisq)
## ulva_hsat_treatment_null
## ulva_hsat_model2 5.948e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ulva_hsat_temperature_null <- lmer(formula = supersat_avg ~ Treatment + (1 | Run) + (1 | Plant.ID) +
(1 | RLC.Order), data = ulva, REML = FALSE)
ulva_hsat_model3 <- lmer(formula = supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +
(1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_hsat_temperature_null, ulva_hsat_model3)
## Data: ulva
## Models:
## ulva_hsat_temperature_null: supersat_avg ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_hsat_model3: supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## npar AIC BIC logLik deviance Chisq Df
## ulva_hsat_temperature_null 9 2393.4 2424.7 -1187.7 2375.4
## ulva_hsat_model3 11 2396.4 2434.7 -1187.2 2374.4 0.9628 2
## Pr(>Chisq)
## ulva_hsat_temperature_null
## ulva_hsat_model3 0.6179
Make plots of the data
ulva %>% ggplot(aes(treatment_graph, supersat_avg)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) +
labs(x="Treatment (salinty/nitrate)", y= "Hsat Time (min)", title= "A", subtitle = "Ulva lactuca") +
scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "18ppt/53umolN", "11ppt/80umolN")) +
ylim(140, 620) + stat_mean() +
geom_hline(yintercept=400, color = "red", size = 0.5, alpha = 0.5) +
scale_color_manual(values = c("#295102", "#7CB950", "#BDE269")) +
theme_bw() +
theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))
ulva %>% ggplot(aes(treatment_graph, deltaEk)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.5, size = 3, aes(color = Temperature), show.legend = FALSE) +
labs(x="Treatment (salinty/nitrate)", y= "delta Ek (%)", title= "A", subtitle = "Ulva lactuca") +
scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "18ppt/53umolN", "11ppt/80umolN")) +
ylim(-80, 120) + stat_mean() +
geom_hline(yintercept=0, color = "orange", size = 0.5, alpha = 0.5) +
scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
theme_bw() +
theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
Summarize the means for variables of interest
ulva %>% group_by(Treatment) %>% summarise_at(vars(supersat_avg), list(mean = mean))
## # A tibble: 5 × 2
## Treatment mean
## <fct> <dbl>
## 1 0 477.
## 2 1 432.
## 3 2 435.
## 4 3 402.
## 5 4 394.
ulva %>% group_by(Run) %>% summarise_at(vars(supersat_avg), list(mean = mean))
## # A tibble: 7 × 2
## Run mean
## <fct> <dbl>
## 1 1 461.
## 2 2 465.
## 3 3 422.
## 4 3.5 352.
## 5 4 349.
## 6 6 452.
## 7 6.5 501.
ulva %>% group_by(Run) %>% summarise_at(vars(day_length_avg), list(mean = mean))
## # A tibble: 7 × 2
## Run mean
## <fct> <dbl>
## 1 1 654.
## 2 2 643.
## 3 3 626.
## 4 3.5 621.
## 5 4 611.
## 6 6 626.
## 7 6.5 633.
ulva %>% group_by(Treatment) %>% summarise_at(vars(deltaEk), list(mean = mean))
## # A tibble: 5 × 2
## Treatment mean
## <fct> <dbl>
## 1 0 -27.8
## 2 1 -22.8
## 3 2 -22.9
## 4 3 1.83
## 5 4 4.61
Add growth rate to this dataset. Open growth dataset, prepare data for insertion
growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/all_runs_growth_011723.csv")
growth_rate$Species <- as.factor(growth_rate$Species)
growth_rate$treatment <- as.factor(growth_rate$treatment)
growth_rate$growth_rate_percent <- (growth_rate$final.weight - growth_rate$Initial.weight) / growth_rate$Initial.weight * 100
gr_ulva <- subset(growth_rate, Species == "Ul" & treatment != 2.5)
ulva$growth_rate <- round((gr_ulva$final.weight - gr_ulva$Initial.weight) / gr_ulva$Initial.weight * 100, digits = 2)
Plot a regression between the photosynthetic independent variables of interest and growth rate
ulva_growth_hsat_graph <- ggplot(ulva, aes(x=supersat_avg, y=growth_rate)) +
geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) +
geom_smooth(method = "lm", col = "black") + theme_bw() +
labs(title = "Ulva lactuca Hsat vs Growth Rate", x = "Hsat time (minutes)",
y = "growth rate (%)") + stat_regline_equation(label.x = 240, label.y = 165) + stat_cor()
ulva_growth_hsat_graph
## `geom_smooth()` using formula = 'y ~ x'
HYPNEA_____________________________________________________________________________________________________________
hypnea <- subset(ek_irrad_data, Species == "hm" & day1_rlc_time != "11:34:05")
hypnea$treatment_graph[hypnea$Treatment == 0] <- "1) 35ppt/0.5umol"
hypnea$treatment_graph[hypnea$Treatment == 1] <- "2) 35ppt/14umol"
hypnea$treatment_graph[hypnea$Treatment == 2] <- "3) 28ppt/27umol"
hypnea$treatment_graph[hypnea$Treatment == 3] <- "5) 18ppt/53umol"
hypnea$treatment_graph[hypnea$Treatment == 4] <- "6) 11ppt/80umol"
hypnea$treatment_graph[hypnea$Treatment == 2.5] <- "4) 28ppt/53umol"
Make a histogram for Hypnea
hist(hypnea$supersat_avg, main = paste("Hypnea musciformis"), col = "maroon", labels = TRUE)
hypnea %>% ggplot(aes(supersat_avg)) +
geom_histogram(binwidth=5, fill = "#a8325e", color = "black", size = 0.25, alpha = 0.85) +
theme_bw()
Run model for hsat hypnea
hsat_model_hypnea <- lmer(formula = supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = hypnea)
Make residual plots of the data
hist(resid(hsat_model_hypnea))
plot(resid(hsat_model_hypnea) ~ fitted(hsat_model_hypnea))
qqnorm(resid(hsat_model_hypnea))
qqline(resid(hsat_model_hypnea))
Check the performance of the model, get r2, plot model and print table of the data
performance::check_model(hsat_model_hypnea)
## Variable `Component` is not in your data frame :/
r.squaredGLMM(hsat_model_hypnea)
## R2m R2c
## [1,] 0.144871 0.7606678
summary(hsat_model_hypnea)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
## Data: hypnea
##
## REML criterion at convergence: 2894.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3056 -0.4446 0.0506 0.5483 2.1875
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant.ID (Intercept) 790.2 28.11
## Run (Intercept) 2049.5 45.27
## Residual 1103.6 33.22
## Number of obs: 286, groups: Plant.ID, 144; Run, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 396.461 32.941 5.042 12.035 6.62e-05 ***
## Treatment1 -49.575 38.936 5.020 -1.273 0.258708
## Treatment2 -35.802 38.936 5.020 -0.920 0.399862
## Treatment2.5 -37.009 56.153 4.730 -0.659 0.540567
## Treatment3 -51.422 38.936 5.020 -1.321 0.243612
## Treatment4 -69.064 38.951 5.027 -1.773 0.136083
## Temperature27 29.379 7.917 131.257 3.711 0.000304 ***
## Temperature30 32.503 7.927 132.169 4.100 7.17e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1 -0.830
## Treatment2 -0.830 0.985
## Treatmnt2.5 -0.575 0.487 0.487
## Treatment3 -0.830 0.985 0.985 0.487
## Treatment4 -0.830 0.984 0.984 0.487 0.984
## Temperatr27 -0.120 0.000 0.000 0.000 0.000 0.000
## Temperatr30 -0.120 0.000 0.000 0.000 0.000 0.001 0.499
tab_model(hsat_model_hypnea, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
| supersat_avg | ||||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | Statistic | p | df |
| (Intercept) | 396.46 | 32.94 | 331.61 – 461.31 | 12.04 | <0.001 | 275.00 |
| Treatment [1] | -49.58 | 38.94 | -126.23 – 27.08 | -1.27 | 0.204 | 275.00 |
| Treatment [2] | -35.80 | 38.94 | -112.45 – 40.85 | -0.92 | 0.359 | 275.00 |
| Treatment [2.5] | -37.01 | 56.15 | -147.55 – 73.53 | -0.66 | 0.510 | 275.00 |
| Treatment [3] | -51.42 | 38.94 | -128.07 – 25.23 | -1.32 | 0.188 | 275.00 |
| Treatment [4] | -69.06 | 38.95 | -145.74 – 7.62 | -1.77 | 0.077 | 275.00 |
| Temperature [27] | 29.38 | 7.92 | 13.79 – 44.97 | 3.71 | <0.001 | 275.00 |
| Temperature [30] | 32.50 | 7.93 | 16.90 – 48.11 | 4.10 | <0.001 | 275.00 |
| Random Effects | ||||||
| σ2 | 1103.65 | |||||
| τ00 Plant.ID | 790.18 | |||||
| τ00 Run | 2049.48 | |||||
| ICC | 0.72 | |||||
| N Run | 8 | |||||
| N Plant.ID | 144 | |||||
| Observations | 286 | |||||
| Marginal R2 / Conditional R2 | 0.145 / 0.761 | |||||
plot(allEffects(hsat_model_hypnea))
Construct null model to perform likelihood ratio test REML must be FALSE
hyp_hsat_treatment_null <- lmer(formula = supersat_avg ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hyp_hsat_model2 <- lmer(formula = supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hyp_hsat_treatment_null, hyp_hsat_model2)
## Data: hypnea
## Models:
## hyp_hsat_treatment_null: supersat_avg ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hyp_hsat_model2: supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## npar AIC BIC logLik deviance Chisq Df
## hyp_hsat_treatment_null 7 2987.0 3012.6 -1486.5 2973.0
## hyp_hsat_model2 12 2972.1 3015.9 -1474.0 2948.1 24.914 5
## Pr(>Chisq)
## hyp_hsat_treatment_null
## hyp_hsat_model2 0.0001448 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hyp_hsat_temperature_null <- lmer(formula = supersat_avg ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hyp_hsat_model3 <- lmer(formula = supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hyp_hsat_temperature_null, hyp_hsat_model3)
## Data: hypnea
## Models:
## hyp_hsat_temperature_null: supersat_avg ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hyp_hsat_model3: supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## npar AIC BIC logLik deviance Chisq Df
## hyp_hsat_temperature_null 10 2985.0 3021.5 -1482.5 2965.0
## hyp_hsat_model3 12 2972.1 3015.9 -1474.0 2948.1 16.894 2
## Pr(>Chisq)
## hyp_hsat_temperature_null
## hyp_hsat_model3 0.0002145 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Make plots
hypnea %>% ggplot(aes(treatment_graph, supersat_avg)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) +
labs(x="Treatment (salinty/nitrate)", y= "Hsat Time (min)", title= "B", subtitle = "Hypnea musciformis") +
scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) +
ylim(140, 620) + stat_mean() +
geom_hline(yintercept=400, color = "red", size = 0.5, alpha = 0.5) +
scale_color_manual(values = c("#9C0627", "#BB589F", "#F4B4E2")) +
theme_bw() +
theme(legend.position = c(0.88,0.88), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))
Plot deltaEk for Hypnea
hypnea %>% ggplot(aes(treatment_graph, deltaEk)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.5, size = 3, aes(color = Temperature), show.legend = FALSE) +
labs(x="Treatment (salinty/nitrate)", y= "delta Ek (%)", title= "B", subtitle = "Hypnea musciformis") +
scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) +
ylim(-100, 250) + stat_mean() +
geom_hline(yintercept=0, color = "orange", size = 0.5, alpha = 0.5) +
scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
theme_bw() +
theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
Prep data for linear regression with growth. For Hypnea remove hm6-4 on 11/12 that had no d9 RLC (final weight 0.1017) and hm6-4 on 10/09/21 because it was white and also looked dead
gr_hypnea <- subset(growth_rate, Species == "Hm" & final.weight != 0.1017 & growth_rate_percent > -87.96837)
hypnea$growth_rate <- round((gr_hypnea$final.weight - gr_hypnea$Initial.weight) / gr_hypnea$Initial.weight * 100, digits = 2)
Plot a regression between the photosynthetic independent variables of interest and growth rate
hypnea_growth_irrad_ek_graph <- ggplot(hypnea, aes(x=supersat_avg, y=growth_rate)) +
geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) +
geom_smooth(method = "lm", col = "black") + theme_bw() +
labs(title = "B", subtitle = "Hypnea musciformis", x = "Hsat Time (min)", y = "growth rate (%)") +
stat_regline_equation(label.x = 400, label.y = 150) + stat_cor(label.x = 400, label.y = 140) +
theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
hypnea_growth_irrad_ek_graph
## `geom_smooth()` using formula = 'y ~ x'
Summarize the means for rETRmax
hypnea %>% group_by(Treatment) %>% summarise_at(vars(supersat_avg), list(mean = mean))
## # A tibble: 6 × 2
## Treatment mean
## <fct> <dbl>
## 1 0 417.
## 2 1 374.
## 3 2 388.
## 4 2.5 380.
## 5 3 373.
## 6 4 354.
hypnea %>% group_by(Run) %>% summarise_at(vars(supersat_avg), list(mean = mean))
## # A tibble: 8 × 2
## Run mean
## <fct> <dbl>
## 1 1 407.
## 2 2 427.
## 3 3 360.
## 4 3.5 312.
## 5 4 320.
## 6 7 380.
## 7 8 406.
## 8 9 428.
hypnea %>% group_by(Run) %>% summarise_at(vars(day_length_avg), list(mean = mean))
## # A tibble: 8 × 2
## Run mean
## <fct> <dbl>
## 1 1 654.
## 2 2 645.
## 3 3 626.
## 4 3.5 621.
## 5 4 611.
## 6 7 688.
## 7 8 635.
## 8 9 636.
hypnea %>% group_by(Treatment) %>% summarise_at(vars(deltaEk), list(mean = mean))
## # A tibble: 6 × 2
## Treatment mean
## <fct> <dbl>
## 1 0 9.46
## 2 1 -37.4
## 3 2 -33.2
## 4 2.5 7.35
## 5 3 -17.1
## 6 4 -21.1